This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.
Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.
# example code from textbook
library(astsa)
par(mfrow=2:1)
tsplot(jj, ylab="QEPS", type="o", col=4, main="Johnson & Johnson
Quarterly Earnings")
tsplot(log(jj), ylab="log(QEPS)", type="o", col=4)
library(xts)
djia_return = diff(log(djia$Close))[-1]
par(mfrow=2:1)
plot(djia$Close, col=4)
plot(djia_return, col=4)
tsplot(diff(log(gdp)), type="o", col=4, ylab="GDP Growth") # diff-log
points(diff(gdp)/lag(gdp,-1), pch=3, col=2)
par(mfrow = c(2,1))
tsplot(soi, ylab="", xlab="", main="Southern Oscillation Index", col=4)
text(1970, .91, "COOL", col="cyan4")
text(1970,-.91, "WARM", col="darkmagenta")
tsplot(rec, ylab="", main="Recruitment", col=4)
culer = c(rgb(.85,.30,.12,.6), rgb(.12,.67,.86,.6))
tsplot(Hare, col = culer[1], lwd=2, type="o", pch=0,
ylab=expression(Number~~~(""%*% 1000)))
lines(Lynx, col=culer[2], lwd=2, type="o", pch=2)
legend("topright", col=culer, lty=1, lwd=2, pch=c(0,2),
legend=c("Hare", "Lynx"), bty="n")
par(mfrow=c(3,1))
culer = c(rgb(.12,.67,.85,.7), rgb(.67,.12,.85,.7))
u = rep(c(rep(.6,16), rep(-.6,16)), 4) # stimulus signal
tsplot(fmri1[,4], ylab="BOLD", xlab="", main="Cortex", col=culer[1],
ylim=c(-.6,.6), lwd=2)
lines(fmri1[,5], col=culer[2], lwd=2)
lines(u, type="s")
tsplot(fmri1[,6], ylab="BOLD", xlab="", main="Thalamus", col=culer[1],
ylim=c(-.6,.6), lwd=2)
lines(fmri1[,7], col=culer[2], lwd=2)
lines(u, type="s")
tsplot(fmri1[,8], ylab="BOLD", xlab="", main="Cerebellum",
col=culer[1], ylim=c(-.6,.6), lwd=2)
lines(fmri1[,9], col=culer[2], lwd=2)
lines(u, type="s")
mtext("Time (1 pt = 2 sec)", side=1, line=1.75)
par(mfrow=2:1)
w = rnorm(500) # 500 N(0,1) variates
v = filter(w, sides=2, filter=rep(1/3,3)) # moving average
tsplot(w, col=4, main="white noise")
tsplot(v, ylim=c(-3,3), col=4, main="moving average")
set.seed(90210)
w = rnorm(250 + 50) # 50 extra to avoid startup problems
x = filter(w, filter=c(1.5,-.75), method="recursive")[-(1:50)]
tsplot(x, main="autoregression", col=4)
set.seed(314159265) # so you can reproduce the results
w = rnorm(200); x = cumsum(w) # random walk
wd = w +.3; xd = cumsum(wd) # random walk with drift
tsplot(xd, ylim=c(-2,80), main="random walk", ylab="", col=4)
abline(a=0, b=.3, lty=2, col=4) # plot drift
lines(x, col="darkred")
abline(h=0, col="darkred", lty=2)
t = 1:500
cs = 2*cos(2*pi*(t+15)/50) # signal
w = rnorm(500) # noise
par(mfrow=c(3,1))
tsplot(cs, col=4, main=expression(2*cos(2*pi*(t+15)/50)))
tsplot(cs+w, col=4, main=expression(2*cos(2*pi*(t+15)/50+N(0,1))))
tsplot(cs+5*w,col=4, main=expression(2*cos(2*pi*(t+15)/50)+N(0,5^2)))
ACF = c(0,0,0,1,2,3,2,1,0,0,0)/3
LAG = -5:5
tsplot(LAG, ACF, type="h", lwd=3, xlab="LAG")
abline(h=0)
points(LAG[-(4:8)], ACF[-(4:8)], pch=20)
axis(1, at=seq(-5, 5, by=2))
x = rnorm(100)
y = lag(x,-5) + rnorm(100)
ccf(y, x, ylab="CCovF", type="covariance", panel.first=Grid())
(r = acf1(soi, 6, plot=FALSE)) # sample acf values
[1] 0.60410089 0.37379533 0.21412447 0.05013659 -0.10703704 -0.18698742
par(mfrow=c(1,2), mar=c(2.5,2.5,0,0)+.5, mgp=c(1.6,.6,0))
plot(lag(soi,-1), soi, col="dodgerblue3", panel.first=Grid())
legend("topleft", legend=r[1], bg="white", adj=.45, cex = 0.85)
plot(lag(soi,-6), soi, col="dodgerblue3", panel.first=Grid())
legend("topleft", legend=r[6], bg="white", adj=.25, cex = 0.8)
set.seed(101011)
x1 = sample(c(-2,2), 11, replace=TRUE) # simulated coin tosses
x2 = sample(c(-2,2), 101, replace=TRUE)
y1 = 5 + filter(x1, sides=1, filter=c(1,-.5))[-1]
y2 = 5 + filter(x2, sides=1, filter=c(1,-.5))[-1]
tsplot(y1, type="s", col=4, xaxt="n", yaxt="n") # y2 not shown
axis(1, 1:10); axis(2, seq(2,8,2), las=1)
points(y1, pch=21, cex=1.1, bg=6)
acf(y1, lag.max=4, plot=FALSE)
Autocorrelations of series ‘y1’, by lag
0 1 2 3 4
1.000 -0.223 -0.623 0.219 0.416
acf(y2, lag.max=4, plot=FALSE)
Autocorrelations of series ‘y2’, by lag
0 1 2 3 4
1.000 -0.435 0.043 0.061 -0.075
par(mfrow=c(3,1))
acf1(soi, 48, main="Southern Oscillation Index")
[1] 0.60 0.37 0.21 0.05 -0.11 -0.19 -0.18 -0.10 0.05 0.22 0.36 0.41 0.31 0.10 -0.06 -0.17 -0.29 -0.37 -0.32 -0.19 -0.04
[22] 0.15 0.31 0.35 0.25 0.10 -0.03 -0.16 -0.28 -0.37 -0.32 -0.16 -0.02 0.17 0.33 0.39 0.30 0.16 0.00 -0.13 -0.24 -0.27
[43] -0.25 -0.13 0.06 0.21 0.38 0.40
acf1(rec, 48, main="Recruitment")
[1] 0.92 0.78 0.63 0.48 0.36 0.26 0.18 0.13 0.09 0.07 0.06 0.02 -0.04 -0.12 -0.19 -0.24 -0.27 -0.27 -0.24 -0.19 -0.11
[22] -0.03 0.03 0.06 0.06 0.02 -0.02 -0.06 -0.09 -0.12 -0.13 -0.11 -0.05 0.02 0.08 0.12 0.10 0.06 0.01 -0.02 -0.03 -0.03
[43] -0.02 0.01 0.06 0.12 0.17 0.20
ccf2(soi, rec, 48, main="SOI & Recruitment")
set.seed(1492)
num = 120
t = 1:num
X = ts( 2*cos(2*pi*t/12) + rnorm(num), freq=12 )
Y = ts( 2*cos(2*pi*(t+5)/12) + rnorm(num), freq=12 )
Yw = resid(lm(Y~ cos(2*pi*t/12) + sin(2*pi*t/12), na.action=NULL))
par(mfrow=c(3,2))
tsplot(X, col=4); tsplot(Y, col=4)
acf1(X, 48); acf1(Y, 48)
[1] 0.51 0.32 -0.03 -0.31 -0.58 -0.59 -0.52 -0.25 0.02 0.36 0.46 0.64 0.51 0.26 -0.09 -0.30 -0.53 -0.58 -0.45 -0.21 0.03
[22] 0.36 0.47 0.53 0.40 0.25 -0.10 -0.30 -0.53 -0.53 -0.41 -0.15 0.07 0.30 0.46 0.51 0.31 0.16 -0.11 -0.31 -0.44 -0.42
[43] -0.36 -0.13 0.12 0.30 0.37 0.41
[1] 0.53 0.33 -0.01 -0.31 -0.47 -0.62 -0.50 -0.34 0.03 0.38 0.53 0.54 0.46 0.25 -0.02 -0.28 -0.51 -0.58 -0.50 -0.27 -0.02
[22] 0.28 0.37 0.49 0.41 0.22 0.03 -0.30 -0.48 -0.52 -0.39 -0.23 -0.02 0.24 0.41 0.48 0.42 0.26 -0.01 -0.18 -0.31 -0.43
[43] -0.38 -0.19 -0.01 0.23 0.34 0.40
ccf2(X, Y, 24); ccf2(X, Yw, 24, ylim=c(-.6,.6))
summary(fit <- lm(salmon~time(salmon), na.action=NULL))
Call:
lm(formula = salmon ~ time(salmon), na.action = NULL)
Residuals:
Min 1Q Median 3Q Max
-1.69187 -0.62453 -0.07024 0.51561 2.34959
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -503.08947 34.44164 -14.61 <2e-16 ***
time(salmon) 0.25290 0.01713 14.76 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.8814 on 164 degrees of freedom
Multiple R-squared: 0.5706, Adjusted R-squared: 0.568
F-statistic: 217.9 on 1 and 164 DF, p-value: < 2.2e-16
tsplot(salmon, col=4, ylab="USD per KG", main="Salmon Export Price")
abline(fit)
culer = c(rgb(.66,.12,.85), rgb(.12,.66,.85), rgb(.85,.30,.12))
par(mfrow=c(3,1))
tsplot(cmort, main="Cardiovascular Mortality", col=culer[1],
type="o", pch=19, ylab="")
tsplot(tempr, main="Temperature", col=culer[2], type="o", pch=19,
ylab="")
tsplot(part, main="Particulates", col=culer[3], type="o", pch=19,
ylab="")
##-- Figure 3.3 --##
tsplot(cmort, main="", ylab="", ylim=c(20,130), col=culer[1])
lines(tempr, col=culer[2])
lines(part, col=culer[3])
legend("topright", legend=c("Mortality", "Temperature", "Pollution"),
lty=1, lwd=2, col=culer, bg="white")
panel.cor <- function(x, y, ...){
usr <- par("usr"); on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r <- round(cor(x, y), 2)
text(0.5, 0.5, r, cex = 1.75)
}
pairs(cbind(Mortality=cmort, Temperature=tempr, Particulates=part),
col="dodgerblue3", lower.panel=panel.cor)
par(mfrow = 2:1)
plot(tempr, tempr^2) # collinear
cor(tempr, tempr^2)
temp = tempr - mean(tempr)
plot(temp, temp^2) # not collinear
cor(temp, temp^2)
temp = tempr - mean(tempr) # center temperature
temp2 = temp^2
trend = time(cmort) # time is trend
fit = lm(cmort~ trend + temp + temp2 + part, na.action=NULL)
summary(fit) # regression results
Call:
lm(formula = cmort ~ trend + temp + temp2 + part, na.action = NULL)
Residuals:
Min 1Q Median 3Q Max
-19.0760 -4.2153 -0.4878 3.7435 29.2448
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.831e+03 1.996e+02 14.19 < 2e-16 ***
trend -1.396e+00 1.010e-01 -13.82 < 2e-16 ***
temp -4.725e-01 3.162e-02 -14.94 < 2e-16 ***
temp2 2.259e-02 2.827e-03 7.99 9.26e-15 ***
part 2.554e-01 1.886e-02 13.54 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 6.385 on 503 degrees of freedom
Multiple R-squared: 0.5954, Adjusted R-squared: 0.5922
F-statistic: 185 on 4 and 503 DF, p-value: < 2.2e-16
summary(aov(fit)) # ANOVA table (compare to next line)
Df Sum Sq Mean Sq F value Pr(>F)
trend 1 10667 10667 261.62 <2e-16 ***
temp 1 8607 8607 211.09 <2e-16 ***
temp2 1 3429 3429 84.09 <2e-16 ***
part 1 7476 7476 183.36 <2e-16 ***
Residuals 503 20508 41
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
summary(aov(lm(cmort~cbind(trend, temp, temp2, part)))) # Table 3.1
Df Sum Sq Mean Sq F value Pr(>F)
cbind(trend, temp, temp2, part) 4 30178 7545 185 <2e-16 ***
Residuals 503 20508 41
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
num = length(cmort) # sample size
AIC(fit)/num - log(2*pi) # AIC
[1] 4.721732
BIC(fit)/num - log(2*pi) # BIC
[1] 4.771699
fish = ts.intersect( rec, soiL6=lag(soi,-6) )
summary(fit1 <- lm(rec~ soiL6, data=fish, na.action=NULL))
Call:
lm(formula = rec ~ soiL6, data = fish, na.action = NULL)
Residuals:
Min 1Q Median 3Q Max
-65.187 -18.234 0.354 16.580 55.790
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 65.790 1.088 60.47 <2e-16 ***
soiL6 -44.283 2.781 -15.92 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 22.5 on 445 degrees of freedom
Multiple R-squared: 0.3629, Adjusted R-squared: 0.3615
F-statistic: 253.5 on 1 and 445 DF, p-value: < 2.2e-16
library(dynlm)
summary(fit2 <- dynlm(rec~ L(soi,6)))
Time series regression with "ts" data:
Start = 1950(7), End = 1987(9)
Call:
dynlm(formula = rec ~ L(soi, 6))
Residuals:
Min 1Q Median 3Q Max
-65.187 -18.234 0.354 16.580 55.790
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 65.790 1.088 60.47 <2e-16 ***
L(soi, 6) -44.283 2.781 -15.92 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 22.5 on 445 degrees of freedom
Multiple R-squared: 0.3629, Adjusted R-squared: 0.3615
F-statistic: 253.5 on 1 and 445 DF, p-value: < 2.2e-16
fit = lm(salmon~time(salmon), na.action=NULL) # the regression
par(mfrow=c(2,1)) # plot transformed data
tsplot(resid(fit), main="detrended salmon price")
tsplot(diff(salmon), main="differenced salmon price")
par(mfrow=c(2,1)) # plot their ACFs
acf1(resid(fit), 48, main="detrended salmon price")
[1] 0.89 0.72 0.56 0.42 0.29 0.20 0.14 0.10 0.08 0.06 0.01 -0.07 -0.20 -0.33 -0.43 -0.50 -0.54 -0.58 -0.56 -0.50 -0.44
[22] -0.35 -0.26 -0.19 -0.17 -0.17 -0.17 -0.17 -0.13 -0.09 0.01 0.14 0.27 0.37 0.45 0.48 0.42 0.34 0.27 0.21 0.17 0.14
[43] 0.14 0.15 0.16 0.14 0.10 0.05
acf1(diff(salmon), 48, main="differenced salmon price")
[1] 0.26 -0.02 -0.08 -0.10 -0.19 -0.11 -0.10 -0.06 0.02 0.07 0.15 0.24 0.00 -0.15 -0.14 -0.07 -0.04 -0.23 -0.15 -0.05 -0.07
[22] -0.06 0.13 0.21 0.11 -0.01 0.00 -0.16 -0.09 -0.23 -0.10 0.04 0.12 0.09 0.23 0.34 0.08 -0.06 -0.06 -0.08 -0.07 -0.09
[43] -0.06 0.06 0.10 0.06 0.04 0.12
par(mfrow=c(2,1))
tsplot(diff(gtemp_land), col=4, main="differenced global temperature")
mean(diff(gtemp_land)) # drift since 1880
[1] 0.01595376
acf1(diff(gtemp_land))
[1] -0.51 0.04 0.00 -0.05 0.00 0.09 -0.06 -0.06 0.13 -0.05 -0.13 0.24 -0.15 0.01 0.07 -0.05 -0.02 0.03 0.05 -0.08 0.03
[22] 0.00 -0.08 0.14
mean(window(diff(gtemp_land), start=1980)) # drift since 1980
[1] 0.04909091
layout(matrix(1:4,2), widths=c(2.5,1))
par(mgp=c(1.6,.6,0), mar=c(2,2,.5,0)+.5)
tsplot(varve, main="", ylab="", col=4, margin=0)
mtext("varve", side=3, line=.5, cex=1.2, font=2, adj=0)
tsplot(log(varve), main="", ylab="", col=4, margin=0)
mtext("log(varve)", side=3, line=.5, cex=1.2, font=2, adj=0)
qqnorm(varve, main="", col=4); qqline(varve, col=2, lwd=2)
qqnorm(log(varve), main="", col=4); qqline(log(varve), col=2, lwd=2)
lag1.plot(soi, 12, col="dodgerblue3") # Figure 3.10
lag2.plot(soi, rec, 8, col="dodgerblue3") # Figure 3.11
library(zoo) # zoo allows easy use of the variable names
dummy = ifelse(soi<0, 0, 1)
fish = as.zoo(ts.intersect(rec, soiL6=lag(soi,-6), dL6=lag(dummy,-6)))
summary(fit <- lm(rec~ soiL6*dL6, data=fish, na.action=NULL))
Call:
lm(formula = rec ~ soiL6 * dL6, data = fish, na.action = NULL)
Residuals:
Min 1Q Median 3Q Max
-63.291 -15.821 2.224 15.791 61.788
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 74.479 2.865 25.998 < 2e-16 ***
soiL6 -15.358 7.401 -2.075 0.0386 *
dL6 -1.139 3.711 -0.307 0.7590
soiL6:dL6 -51.244 9.523 -5.381 1.2e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 21.84 on 443 degrees of freedom
Multiple R-squared: 0.4024, Adjusted R-squared: 0.3984
F-statistic: 99.43 on 3 and 443 DF, p-value: < 2.2e-16
plot(fish$soiL6, fish$rec, panel.first=Grid(), col="dodgerblue3")
points(fish$soiL6, fitted(fit), pch=3, col=6)
lines(lowess(fish$soiL6, fish$rec), col=4, lwd=2)
tsplot(resid(fit)) # not shown, but looks like Figure 3.5
acf1(resid(fit)) # and obviously not noise
[1] 0.69 0.62 0.49 0.37 0.24 0.15 0.08 0.00 -0.03 -0.10 -0.13 -0.16 -0.17 -0.23 -0.24 -0.23 -0.23 -0.22 -0.17 -0.09 -0.05
[22] 0.01 0.05 0.06 0.09 0.07 0.10 0.06 0.02 -0.02 -0.02 -0.02
set.seed(90210) # so you can reproduce these results
x = 2*cos(2*pi*1:500/50 + .6*pi) + rnorm(500,0,5)
z1 = cos(2*pi*1:500/50)
z2 = sin(2*pi*1:500/50)
summary(fit <- lm(x~ 0 + z1 + z2)) # zero to exclude the intercept
Call:
lm(formula = x ~ 0 + z1 + z2)
Residuals:
Min 1Q Median 3Q Max
-14.8584 -3.8525 -0.3186 3.3487 15.5440
Coefficients:
Estimate Std. Error t value Pr(>|t|)
z1 -0.7442 0.3274 -2.273 0.0235 *
z2 -1.9949 0.3274 -6.093 2.23e-09 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 5.177 on 498 degrees of freedom
Multiple R-squared: 0.07827, Adjusted R-squared: 0.07456
F-statistic: 21.14 on 2 and 498 DF, p-value: 1.538e-09
par(mfrow=c(2,1))
tsplot(x, col=4)
tsplot(x, ylab=expression(hat(x)), col=rgb(0,0,1,.5))
lines(fitted(fit), col=2, lwd=2)
w = c(.5, rep(1,11), .5)/12
soif = filter(soi, sides=2, filter=w)
tsplot(soi, col=rgb(.5, .6, .85, .9), ylim=c(-1, 1.15))
lines(soif, lwd=2, col=4)
# insert
par(fig = c(.65, 1, .75, 1), new = TRUE)
w1 = c(rep(0,20), w, rep(0,20))
plot(w1, type="l", ylim = c(-.02,.1), xaxt="n", yaxt="n", ann=FALSE)
tsplot(soi, col=rgb(0.5, 0.6, 0.85, .9), ylim=c(-1, 1.15))
lines(ksmooth(time(soi), soi, "normal", bandwidth=1), lwd=2, col=4)
# insert
par(fig = c(.65, 1, .75, 1), new = TRUE)
curve(dnorm(x), -3, 3, xaxt="n", yaxt="n", ann=FALSE, col=4)
SOI = ts(soi, freq=1)
tsplot(SOI) # the time scale matters (not shown)
lines(ksmooth(time(SOI), SOI, "normal", bandwidth=12), lwd=2, col=4)
tsplot(soi, col=rgb(0.5, 0.6, 0.85, .9))
lines(lowess(soi, f=.05), lwd=2, col=4) # El Niño cycle
# lines(lowess(soi), lty=2, lwd=2, col=2) # trend (with default span)
##-- trend with CIs using loess --##
lo = predict(loess(soi~ time(soi)), se=TRUE)
trnd = ts(lo$fit, start=1950, freq=12) # put back ts attributes
lines(trnd, col=6, lwd=2)
L = trnd - qt(0.975, lo$df)*lo$se
U = trnd + qt(0.975, lo$df)*lo$se
xx = c(time(soi), rev(time(soi)))
yy = c(L, rev(U))
polygon(xx, yy, border=8, col=gray(.6, alpha=.4) )
plot(tempr, cmort, xlab="Temperature", ylab="Mortality",
col="dodgerblue3", panel.first=Grid())
lines(lowess(tempr,cmort), col=4, lwd=2)
x = window(hor, start=2002)
plot(decompose(x)) # not shown
plot(stl(x, s.window="per")) # seasons are periodic - not shown
plot(stl(x, s.window=15))
culer = c("cyan4", 4, 2, 6)
par(mfrow = c(4,1), cex.main=1)
x = window(hor, start=2002)
out = stl(x, s.window=15)$time.series
tsplot(x, main="Hawaiian Occupancy Rate", ylab="% rooms", col=gray(.7))
text(x, labels=1:4, col=culer, cex=1.25)
tsplot(out[,1], main="Seasonal", ylab="% rooms",col=gray(.7))
text(out[,1], labels=1:4, col=culer, cex=1.25)
tsplot(out[,2], main="Trend", ylab="% rooms", col=gray(.7))
text(out[,2], labels=1:4, col=culer, cex=1.25)
tsplot(out[,3], main="Noise", ylab="% rooms", col=gray(.7))
text(out[,3], labels=1:4, col=culer, cex=1.25)
par(mfrow=c(2,1))
tsplot(arima.sim(list(order=c(1,0,0), ar=.9), n=100), ylab="x", col=4,
main=expression(AR(1)~~~phi==+.9))
tsplot(arima.sim(list(order=c(1,0,0),ar=-.9), n=100), ylab="x", col=4,
main=expression(AR(1)~~~phi==-.9))
psi = ARMAtoMA(ar = c(1.5, -.75), ma = 0, 50)
par(mfrow=c(2,1), mar=c(2,2.5,1,0)+.5, mgp=c(1.5,.6,0), cex.main=1.1)
plot(psi, xaxp=c(0,144,12), type="n", col=4,
ylab=expression(psi-weights),
main=expression(AR(2)~~~phi[1]==1.5~~~phi[2]==-.75))
abline(v=seq(0,48,by=12), h=seq(-.5,1.5,.5), col=gray(.9))
lines(psi, type="o", col=4)
set.seed(8675309)
simulation = arima.sim(list(order=c(2,0,0), ar=c(1.5,-.75)), n=144)
plot(simulation, xaxp=c(0,144,12), type="n", ylab=expression(X[~t]))
abline(v=seq(0,144,by=12), h=c(-5,0,5), col=gray(.9))
lines(simulation, col=4)
par(mfrow = c(2,1))
tsplot(arima.sim(list(order=c(0,0,1), ma=.9), n=100), col=4,
ylab="x", main=expression(MA(1)~~~theta==+.5))
tsplot(arima.sim(list(order=c(0,0,1), ma=-.9), n=100), col=4,
ylab="x", main=expression(MA(1)~~~theta==-.5))
set.seed(8675309) # Jenny, I got your number
x = rnorm(150, mean=5) # generate iid N(5,1)s
arima(x, order=c(1,0,1)) # estimation
Call:
arima(x = x, order = c(1, 0, 1))
Coefficients:
ar1 ma1 intercept
-0.9595 0.9527 5.0462
s.e. 0.1688 0.1750 0.0727
sigma^2 estimated as 0.7986: log likelihood = -195.98, aic = 399.96
AR = c(1, -.3, -.4) # original AR coefs on the left
polyroot(AR)
[1] 1.25-0i -2.00+0i
MA = c(1, .5) # original MA coefs on the right
polyroot(MA)
[1] -2+0i
ACF = ARMAacf(ar=c(1.5,-.75), ma=0, 50)
plot(ACF, type="h", xlab="lag", panel.first=Grid())
abline(h=0)
ACF = ARMAacf(ar=c(1.5,-.75), ma=0, 24)[-1]
PACF = ARMAacf(ar=c(1.5,-.75), ma=0, 24, pacf=TRUE)
par(mfrow=1:2)
tsplot(ACF, type="h", xlab="lag", ylim=c(-.8,1))
abline(h=0)
tsplot(PACF, type="h", xlab="lag", ylim=c(-.8,1))
abline(h=0)
acf2(rec, 48) # will produce values and a graphic
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21]
ACF 0.92 0.78 0.63 0.48 0.36 0.26 0.18 0.13 0.09 0.07 0.06 0.02 -0.04 -0.12 -0.19 -0.24 -0.27 -0.27 -0.24 -0.19 -0.11
PACF 0.92 -0.44 -0.05 -0.02 0.07 -0.03 -0.03 0.04 0.05 -0.02 -0.05 -0.14 -0.15 -0.05 0.05 0.01 0.01 0.02 0.09 0.11 0.03
[,22] [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37] [,38] [,39] [,40] [,41]
ACF -0.03 0.03 0.06 0.06 0.02 -0.02 -0.06 -0.09 -0.12 -0.13 -0.11 -0.05 0.02 0.08 0.12 0.10 0.06 0.01 -0.02 -0.03
PACF -0.03 -0.01 -0.07 -0.12 -0.03 0.05 -0.08 -0.04 -0.03 0.06 0.05 0.15 0.09 -0.04 -0.10 -0.09 -0.02 0.05 0.08 -0.02
[,42] [,43] [,44] [,45] [,46] [,47] [,48]
ACF -0.03 -0.02 0.01 0.06 0.12 0.17 0.20
PACF -0.01 -0.02 0.05 0.01 0.05 0.08 -0.04
(regr = ar.ols(rec, order=2, demean=FALSE, intercept=TRUE))
Call:
ar.ols(x = rec, order.max = 2, demean = FALSE, intercept = TRUE)
Coefficients:
1 2
1.3541 -0.4632
Intercept: 6.737 (1.111)
Order selected 2 sigma^2 estimated as 89.72
regr$asy.se.coef # standard errors of the estimates
$x.mean
[1] 1.110599
$ar
[1] 0.04178901 0.04187942
set.seed(2)
ma1 = arima.sim(list(order = c(0,0,1), ma = 0.9), n = 50)
acf1(ma1, plot=FALSE)[1]
[1] 0.5066599
tsplot(diff(log(varve)), col=4, ylab=expression(nabla~log~X[~t]),
main="Transformed Glacial Varves")
acf2(diff(log(varve)))
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21]
ACF -0.4 -0.04 -0.06 0.01 0.00 0.04 -0.04 0.04 0.01 -0.05 0.06 -0.06 -0.04 0.08 -0.02 0.01 0.00 0.03 -0.05 -0.06 0.07
PACF -0.4 -0.24 -0.23 -0.18 -0.15 -0.08 -0.11 -0.05 -0.01 -0.07 0.02 -0.05 -0.12 -0.03 -0.05 -0.04 -0.03 0.03 -0.03 -0.13 -0.04
[,22] [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36]
ACF 0.04 -0.06 0.05 -0.01 -0.04 0.05 -0.05 0.03 -0.02 0.00 0.06 -0.05 -0.03 0.04 -0.05
PACF 0.01 -0.06 0.01 0.02 -0.04 0.03 -0.02 0.00 -0.03 -0.02 0.04 -0.02 -0.02 0.01 -0.07
x = diff(log(varve)) # data
r = acf1(x, 1, plot=FALSE) # acf(1)
c(0) -> w -> z -> Sc -> Sz -> Szw -> para # initialize
num = length(x) # = 633
## Estimation
para[1] = (1-sqrt(1-4*(r^2)))/(2*r) # MME
niter = 12
for (j in 1:niter){
for (i in 2:num){ w[i] = x[i] - para[j]*w[i-1]
z[i] = w[i-1] - para[j]*z[i-1]
}
Sc[j] = sum(w^2)
Sz[j] = sum(z^2)
Szw[j] = sum(z*w)
para[j+1] = para[j] + Szw[j]/Sz[j]
}
## Results
cbind(iteration=1:niter-1, thetahat=para[1:niter], Sc, Sz)
iteration thetahat Sc Sz
[1,] 0 -0.4946886 158.7393 171.2397
[2,] 1 -0.6681759 150.7468 235.2660
[3,] 2 -0.7333163 149.2644 300.5618
[4,] 3 -0.7563893 149.0309 336.8234
[5,] 4 -0.7655639 148.9897 354.1729
[6,] 5 -0.7694700 148.9818 362.1665
[7,] 6 -0.7711856 148.9803 365.8014
[8,] 7 -0.7719497 148.9800 367.4456
[9,] 8 -0.7722921 148.9799 368.1875
[10,] 9 -0.7724460 148.9799 368.5220
[11,] 10 -0.7725153 148.9799 368.6727
[12,] 11 -0.7725465 148.9799 368.7406
## Plot conditional SS
c(0) -> w -> cSS
th = -seq(.3, .94, .01)
for (p in 1:length(th)){
for (i in 2:num){ w[i] = x[i] - th[p]*w[i-1]
}
cSS[p] = sum(w^2)
}
tsplot(th, cSS, ylab=expression(S[c](theta)), xlab=expression(theta))
abline(v=para[1:12], lty=2, col=4) # add previous results to plot
points(para[1:12], Sc[1:12], pch=16, col=4)
sarima(diff(log(varve)), p=0, d=0, q=1, no.constant=TRUE)
initial value -0.551778
iter 2 value -0.671626
iter 3 value -0.705973
iter 4 value -0.707314
iter 5 value -0.722372
iter 6 value -0.722738
iter 7 value -0.723187
iter 8 value -0.723194
iter 9 value -0.723195
iter 9 value -0.723195
iter 9 value -0.723195
final value -0.723195
converged
initial value -0.722700
iter 2 value -0.722702
iter 3 value -0.722702
iter 3 value -0.722702
iter 3 value -0.722702
final value -0.722702
converged
<><><><><><><><><><><><><><>
Coefficients:
Estimate SE t.value p.value
ma1 -0.7705 0.0341 -22.6161 0
sigma^2 estimated as 0.2353156 on 632 degrees of freedom
AIC = 1.398791 AICc = 1.398802 BIC = 1.412853
sarima(rec, p=2, d=0, q=0)
initial value 3.332380
iter 2 value 3.251366
iter 3 value 2.564654
iter 4 value 2.430141
iter 5 value 2.258212
iter 6 value 2.253343
iter 7 value 2.248346
iter 8 value 2.248345
iter 9 value 2.248345
iter 10 value 2.248341
iter 11 value 2.248332
iter 12 value 2.248331
iter 13 value 2.248330
iter 13 value 2.248330
iter 13 value 2.248330
final value 2.248330
converged
initial value 2.248862
iter 2 value 2.248857
iter 3 value 2.248855
iter 4 value 2.248855
iter 5 value 2.248854
iter 6 value 2.248854
iter 7 value 2.248854
iter 8 value 2.248854
iter 9 value 2.248853
iter 10 value 2.248853
iter 10 value 2.248853
iter 10 value 2.248853
final value 2.248853
converged
<><><><><><><><><><><><><><>
Coefficients:
Estimate SE t.value p.value
ar1 1.3512 0.0416 32.4933 0
ar2 -0.4612 0.0417 -11.0687 0
xmean 61.8585 4.0039 15.4494 0
sigma^2 estimated as 89.33436 on 450 degrees of freedom
AIC = 7.353244 AICc = 7.353362 BIC = 7.389587
sarima.for(rec, n.ahead=24, p=2, d=0, q=0)
$pred
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
1987 20.36547 26.08036 32.65148
1988 38.89474 44.30006 48.72437 52.20958 54.87831 56.87693 58.34666 59.41079 60.17081 60.70697 61.08091 61.33890
1989 61.51504 61.63405 61.71363 61.76626 61.80068 61.82291 61.83707 61.84596 61.85143
$se
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
1987 9.451686 15.888378 20.464325
1988 23.492457 25.393693 26.537088 27.199368 27.570234 27.771616 27.877923 27.932597 27.960042 27.973508 27.979974 27.983014
1989 27.984413 27.985045 27.985323 27.985444 27.985495 27.985515 27.985524 27.985527 27.985528
abline(h=61.8585, col=4) # display estimated mean
sarima(diff(log(varve)), p=0, d=0, q=1, no.constant=TRUE)
initial value -0.551778
iter 2 value -0.671626
iter 3 value -0.705973
iter 4 value -0.707314
iter 5 value -0.722372
iter 6 value -0.722738
iter 7 value -0.723187
iter 8 value -0.723194
iter 9 value -0.723195
iter 9 value -0.723195
iter 9 value -0.723195
final value -0.723195
converged
initial value -0.722700
iter 2 value -0.722702
iter 3 value -0.722702
iter 3 value -0.722702
iter 3 value -0.722702
final value -0.722702
converged
<><><><><><><><><><><><><><>
Coefficients:
Estimate SE t.value p.value
ma1 -0.7705 0.0341 -22.6161 0
sigma^2 estimated as 0.2353156 on 632 degrees of freedom
AIC = 1.398791 AICc = 1.398802 BIC = 1.412853
sarima(log(varve), p=0, d=1, q=1, no.constant=TRUE)
initial value -0.551778
iter 2 value -0.671626
iter 3 value -0.705973
iter 4 value -0.707314
iter 5 value -0.722372
iter 6 value -0.722738
iter 7 value -0.723187
iter 8 value -0.723194
iter 9 value -0.723195
iter 9 value -0.723195
iter 9 value -0.723195
final value -0.723195
converged
initial value -0.722700
iter 2 value -0.722702
iter 3 value -0.722702
iter 3 value -0.722702
iter 3 value -0.722702
final value -0.722702
converged
<><><><><><><><><><><><><><>
Coefficients:
Estimate SE t.value p.value
ma1 -0.7705 0.0341 -22.6161 0
sigma^2 estimated as 0.2353156 on 632 degrees of freedom
AIC = 1.398792 AICc = 1.398802 BIC = 1.412853
ARMAtoMA(ar=1, ma=0, 20) # ψ-weights
[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
round( ARMAtoMA(ar=c(1.9,-.9), ma=0, 60), 1 )
[1] 1.9 2.7 3.4 4.1 4.7 5.2 5.7 6.1 6.5 6.9 7.2 7.5 7.7 7.9 8.1 8.3 8.5 8.6 8.8 8.9 9.0 9.1 9.2 9.3 9.4
[26] 9.4 9.5 9.5 9.6 9.6 9.7 9.7 9.7 9.7 9.8 9.8 9.8 9.8 9.9 9.9 9.9 9.9 9.9 9.9 9.9 9.9 9.9 9.9 9.9 10.0
[51] 10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0
set.seed(1998)
x <- ts(arima.sim(list(order = c(1,1,0), ar=.9), n=150)[-1])
y <- window(x, start=1, end=100)
sarima.for(y, n.ahead = 50, p = 1, d = 1, q = 0, plot.all=TRUE)
$pred
Time Series:
Start = 101
End = 150
Frequency = 1
[1] 108.8050 107.6806 106.8693 106.3241 106.0053 105.8791 105.9165 106.0932 106.3883 106.7841 107.2656 107.8199 108.4361 109.1050
[15] 109.8187 110.5705 111.3547 112.1665 113.0017 113.8568 114.7289 115.6154 116.5141 117.4233 118.3413 119.2669 120.1989 121.1363
[29] 122.0784 123.0244 123.9738 124.9260 125.8807 126.8374 127.7959 128.7558 129.7171 130.6794 131.6426 132.6066 133.5713 134.5365
[43] 135.5022 136.4684 137.4348 138.4016 139.3686 140.3358 141.3033 142.2708
$se
Time Series:
Start = 101
End = 150
Frequency = 1
[1] 0.9570915 2.0131097 3.1812332 4.4084707 5.6617584 6.9197429 8.1683744 9.3983950 10.6037992 11.7808348 12.9273310
[12] 14.0422326 15.1252705 16.1767257 17.1972588 18.1877853 19.1493851 20.0832363 20.9905665 21.8726189 22.7306265 23.5657957
[23] 24.3792934 25.1722400 25.9457044 26.7007014 27.4381914 28.1590801 28.8642201 29.5544126 30.2304096 30.8929162 31.5425932
[34] 32.1800597 32.8058954 33.4206433 34.0248122 34.6188788 35.2032899 35.7784645 36.3447959 36.9026532 37.4523833 37.9943122
[45] 38.5287468 39.0559758 39.5762716 40.0898907 40.5970755 41.0980549
text(85, 205, "PAST"); text(115, 205, "FUTURE")
abline(v=100, lty=2, col=4)
lines(x)
set.seed(666)
x = arima.sim(list(order = c(0,1,1), ma = -0.8), n = 100)
(x.ima = HoltWinters(x, beta=FALSE, gamma=FALSE)) # α below is 1 − λ
Holt-Winters exponential smoothing without trend and without seasonal component.
Call:
HoltWinters(x = x, beta = FALSE, gamma = FALSE)
Smoothing parameters:
alpha: 0.1663072
beta : FALSE
gamma: FALSE
Coefficients:
[,1]
a -2.241533
plot(x.ima, main="EWMA")
##-- Figure 5.3 --##
layout(1:2, heights=2:1)
tsplot(gnp, col=4)
acf1(gnp, main="")
[1] 0.99 0.97 0.96 0.94 0.93 0.91 0.90 0.88 0.87 0.85 0.83 0.82 0.80 0.79 0.77 0.76 0.74 0.73 0.72 0.70 0.69 0.68 0.66 0.65 0.64
##-- Figure 5.4 --##
tsplot(diff(log(gnp)), ylab="GNP Growth Rate", col=4)
sarima(diff(log(gnp)), 0,0,2) # MA(2) on growth rate
initial value -4.591629
iter 2 value -4.661095
iter 3 value -4.662220
iter 4 value -4.662243
iter 5 value -4.662243
iter 6 value -4.662243
iter 6 value -4.662243
iter 6 value -4.662243
final value -4.662243
converged
initial value -4.662022
iter 2 value -4.662023
iter 2 value -4.662023
iter 2 value -4.662023
final value -4.662023
converged
<><><><><><><><><><><><><><>
Coefficients:
Estimate SE t.value p.value
ma1 0.3028 0.0654 4.6272 0.0000
ma2 0.2035 0.0644 3.1594 0.0018
xmean 0.0083 0.0010 8.7178 0.0000
sigma^2 estimated as 8.919178e-05 on 219 degrees of freedom
AIC = -6.450133 AICc = -6.449637 BIC = -6.388823
sarima(diff(log(gnp)), 1,0,0)
initial value -4.589567
iter 2 value -4.654150
iter 3 value -4.654150
iter 4 value -4.654151
iter 4 value -4.654151
iter 4 value -4.654151
final value -4.654151
converged
initial value -4.655919
iter 2 value -4.655921
iter 3 value -4.655922
iter 4 value -4.655922
iter 5 value -4.655922
iter 5 value -4.655922
iter 5 value -4.655922
final value -4.655922
converged
<><><><><><><><><><><><><><>
Coefficients:
Estimate SE t.value p.value
ar1 0.3467 0.0627 5.5255 0
xmean 0.0083 0.0010 8.5398 0
sigma^2 estimated as 9.029569e-05 on 220 degrees of freedom
AIC = -6.44694 AICc = -6.446693 BIC = -6.400958
round( ARMAtoMA(ar=.35, ma=0, 10), 3) # print psi-weights
[1] 0.350 0.122 0.043 0.015 0.005 0.002 0.001 0.000 0.000 0.000
sarima(diff(log(gnp)), 0, 0, 2) # MA(2) fit with diagnostics
initial value -4.591629
iter 2 value -4.661095
iter 3 value -4.662220
iter 4 value -4.662243
iter 5 value -4.662243
iter 6 value -4.662243
iter 6 value -4.662243
iter 6 value -4.662243
final value -4.662243
converged
initial value -4.662022
iter 2 value -4.662023
iter 2 value -4.662023
iter 2 value -4.662023
final value -4.662023
converged
<><><><><><><><><><><><><><>
Coefficients:
Estimate SE t.value p.value
ma1 0.3028 0.0654 4.6272 0.0000
ma2 0.2035 0.0644 3.1594 0.0018
xmean 0.0083 0.0010 8.7178 0.0000
sigma^2 estimated as 8.919178e-05 on 219 degrees of freedom
AIC = -6.450133 AICc = -6.449637 BIC = -6.388823
sarima(diff(log(gnp)), 0, 0, 3) # try an MA(2+1) fit (not shown)
initial value -4.591629
iter 2 value -4.660639
iter 3 value -4.665404
iter 4 value -4.665858
iter 5 value -4.665953
iter 6 value -4.665955
iter 7 value -4.665956
iter 8 value -4.665956
iter 8 value -4.665956
iter 8 value -4.665956
final value -4.665956
converged
initial value -4.665676
iter 2 value -4.665678
iter 3 value -4.665678
iter 3 value -4.665678
iter 3 value -4.665678
final value -4.665678
converged
<><><><><><><><><><><><><><>
Coefficients:
Estimate SE t.value p.value
ma1 0.3208 0.0662 4.8430 0.0000
ma2 0.2478 0.0718 3.4512 0.0007
ma3 0.0909 0.0701 1.2962 0.1963
xmean 0.0083 0.0010 7.9444 0.0000
sigma^2 estimated as 8.852732e-05 on 218 degrees of freedom
AIC = -6.448434 AICc = -6.447604 BIC = -6.371797
sarima(diff(log(gnp)), 2, 0, 0) # try an AR(1+1) fit (not shown)
initial value -4.588072
iter 2 value -4.647096
iter 3 value -4.655676
iter 4 value -4.655896
iter 5 value -4.655903
iter 6 value -4.655903
iter 7 value -4.655904
iter 7 value -4.655904
iter 7 value -4.655904
final value -4.655904
converged
initial value -4.659296
iter 2 value -4.659300
iter 3 value -4.659306
iter 4 value -4.659307
iter 5 value -4.659307
iter 5 value -4.659307
iter 5 value -4.659307
final value -4.659307
converged
<><><><><><><><><><><><><><>
Coefficients:
Estimate SE t.value p.value
ar1 0.3180 0.0667 4.7660 0.0000
ar2 0.0820 0.0668 1.2282 0.2207
xmean 0.0083 0.0011 7.8742 0.0000
sigma^2 estimated as 8.968109e-05 on 219 degrees of freedom
AIC = -6.444701 AICc = -6.444205 BIC = -6.383392
sarima(log(varve), 0, 1, 1, no.constant=TRUE) # ARIMA(0,1,1)
initial value -0.551778
iter 2 value -0.671626
iter 3 value -0.705973
iter 4 value -0.707314
iter 5 value -0.722372
iter 6 value -0.722738
iter 7 value -0.723187
iter 8 value -0.723194
iter 9 value -0.723195
iter 9 value -0.723195
iter 9 value -0.723195
final value -0.723195
converged
initial value -0.722700
iter 2 value -0.722702
iter 3 value -0.722702
iter 3 value -0.722702
iter 3 value -0.722702
final value -0.722702
converged
<><><><><><><><><><><><><><>
Coefficients:
Estimate SE t.value p.value
ma1 -0.7705 0.0341 -22.6161 0
sigma^2 estimated as 0.2353156 on 632 degrees of freedom
AIC = 1.398792 AICc = 1.398802 BIC = 1.412853
sarima(log(varve), 1, 1, 1, no.constant=TRUE) # ARIMA(1,1,1)
initial value -0.550992
iter 2 value -0.648952
iter 3 value -0.676952
iter 4 value -0.699136
iter 5 value -0.724481
iter 6 value -0.726964
iter 7 value -0.734257
iter 8 value -0.735999
iter 9 value -0.737045
iter 10 value -0.737381
iter 11 value -0.737469
iter 12 value -0.737473
iter 13 value -0.737473
iter 14 value -0.737473
iter 14 value -0.737473
iter 14 value -0.737473
final value -0.737473
converged
initial value -0.737355
iter 2 value -0.737361
iter 3 value -0.737362
iter 4 value -0.737363
iter 5 value -0.737363
iter 5 value -0.737363
iter 5 value -0.737363
final value -0.737363
converged
<><><><><><><><><><><><><><>
Coefficients:
Estimate SE t.value p.value
ar1 0.2330 0.0518 4.4994 0
ma1 -0.8858 0.0292 -30.3861 0
sigma^2 estimated as 0.2284339 on 631 degrees of freedom
AIC = 1.37263 AICc = 1.372661 BIC = 1.393723
uspop = c(75.995, 91.972, 105.711, 123.203, 131.669,150.697,
179.323, 203.212, 226.505, 249.633, 281.422, 308.745)
uspop = ts(uspop, start=1900, freq=.1)
t = time(uspop) - 1955
reg = lm( uspop~ t+I(t^2)+I(t^3)+I(t^4)+I(t^5)+I(t^6)+I(t^7)+I(t^8) )
b = as.vector(reg$coef)
g = function(t){ b[1] + b[2]*(t-1955) + b[3]*(t-1955)^2 +
b[4]*(t-1955)^3 + b[5]*(t-1955)^4 + b[6]*(t-1955)^5 +
b[7]*(t-1955)^6 + b[8]*(t-1955)^7 + b[9]*(t-1955)^8
}
par(mar=c(2,2.5,.5,0)+.5, mgp=c(1.6,.6,0))
curve(g, 1900, 2024, ylab="Population", xlab="Year", main="U.S.
Population by Official Census", panel.first=Grid(),
cex.main=1, font.main=1, col=4)
abline(v=seq(1910,2020,by=20), lty=1, col=gray(.9))
points(time(uspop), uspop, pch=21, bg=rainbow(12), cex=1.25)
mtext(expression(""%*% 10^6), side=2, line=1.5, adj=.95)
axis(1, seq(1910,2020,by=20), labels=TRUE)
sarima(diff(log(gnp)), 1, 0, 0) # AR(1)
initial value -4.589567
iter 2 value -4.654150
iter 3 value -4.654150
iter 4 value -4.654151
iter 4 value -4.654151
iter 4 value -4.654151
final value -4.654151
converged
initial value -4.655919
iter 2 value -4.655921
iter 3 value -4.655922
iter 4 value -4.655922
iter 5 value -4.655922
iter 5 value -4.655922
iter 5 value -4.655922
final value -4.655922
converged
<><><><><><><><><><><><><><>
Coefficients:
Estimate SE t.value p.value
ar1 0.3467 0.0627 5.5255 0
xmean 0.0083 0.0010 8.5398 0
sigma^2 estimated as 9.029569e-05 on 220 degrees of freedom
AIC = -6.44694 AICc = -6.446693 BIC = -6.400958
sarima(diff(log(gnp)), 0, 0, 2) # MA(2)
initial value -4.591629
iter 2 value -4.661095
iter 3 value -4.662220
iter 4 value -4.662243
iter 5 value -4.662243
iter 6 value -4.662243
iter 6 value -4.662243
iter 6 value -4.662243
final value -4.662243
converged
initial value -4.662022
iter 2 value -4.662023
iter 2 value -4.662023
iter 2 value -4.662023
final value -4.662023
converged
<><><><><><><><><><><><><><>
Coefficients:
Estimate SE t.value p.value
ma1 0.3028 0.0654 4.6272 0.0000
ma2 0.2035 0.0644 3.1594 0.0018
xmean 0.0083 0.0010 8.7178 0.0000
sigma^2 estimated as 8.919178e-05 on 219 degrees of freedom
AIC = -6.450133 AICc = -6.449637 BIC = -6.388823
set.seed(666)
phi = c(rep(0,11),.9)
sAR = ts(arima.sim(list(order=c(12,0,0), ar=phi), n=37), freq=12) + 50
layout(matrix(c(1,2, 1,3), nc=2), heights=c(1.5,1))
par(mar=c(2.5,2.5,2,1), mgp=c(1.6,.6,0))
plot(sAR, xaxt="n", col=gray(.6), main="seasonal AR(1)", xlab="YEAR",
type="c", ylim=c(45,54))
abline(v=1:4, lty=2, col=gray(.6))
axis(1,1:4); box()
abline(h=seq(46,54,by=2), col=gray(.9))
Months = c("J","F","M","A","M","J","J","A","S","O","N","D")
points(sAR, pch=Months, cex=1.35, font=4, col=1:4)
ACF = ARMAacf(ar=phi, ma=0, 100)[-1]
PACF = ARMAacf(ar=phi, ma=0, 100, pacf=TRUE)
LAG = 1:100/12
plot(LAG, ACF, type="h", xlab="LAG", ylim=c(-.1,1), axes=FALSE)
segments(0,0,0,1)
axis(1, seq(0,8,by=1)); axis(2); box(); abline(h=0)
plot(LAG, PACF, type="h", xlab="LAG", ylim=c(-.1,1), axes=FALSE)
axis(1, seq(0,8,by=1)); axis(2); box(); abline(h=0)
##-- Figure 5.10 --##
phi = c(rep(0,11),.8)
ACF = ARMAacf(ar=phi, ma=-.5, 50)[-1]
PACF = ARMAacf(ar=phi, ma=-.5, 50, pacf=TRUE)
LAG = 1:50/12
par(mfrow=c(1,2))
plot(LAG, ACF, type="h", ylim=c(-.4,.8), panel.first=Grid())
abline(h=0)
plot(LAG, PACF, type="h", ylim=c(-.4,.8), panel.first=Grid())
abline(h=0)
##-- birth series --##
tsplot(birth) # monthly number of births in US
acf2( diff(birth) ) # P/ACF of the differenced birth rate
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21]
ACF -0.32 0.16 -0.08 -0.19 0.09 -0.28 0.06 -0.19 -0.05 0.17 -0.26 0.82 -0.28 0.17 -0.07 -0.18 0.08 -0.28 0.07 -0.18 -0.05
PACF -0.32 0.06 -0.01 -0.25 -0.03 -0.26 -0.17 -0.29 -0.35 -0.16 -0.59 0.57 0.13 0.11 0.13 0.09 0.00 0.00 0.05 0.04 -0.07
[,22] [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37] [,38] [,39] [,40] [,41]
ACF 0.16 -0.24 0.78 -0.27 0.19 -0.08 -0.17 0.07 -0.29 0.07 -0.15 -0.04 0.14 -0.24 0.75 -0.23 0.16 -0.08 -0.15 0.05
PACF -0.10 -0.20 0.19 0.01 0.05 0.07 0.07 -0.02 -0.06 -0.02 0.09 0.03 -0.06 -0.16 0.03 0.08 -0.10 -0.03 0.07 -0.04
[,42] [,43] [,44] [,45] [,46] [,47] [,48]
ACF -0.25 0.06 -0.18 -0.03 0.15 -0.22 0.72
PACF 0.06 0.04 -0.07 -0.06 0.02 -0.04 0.10
x = window(hor, start=2002)
par(mfrow = c(2,1))
tsplot(x, main="Hawaiian Quarterly Occupancy Rate", ylab=" % rooms",
ylim=c(62,86), col=gray(.7))
text(x, labels=1:4, col=c(3,4,2,6), cex=.8)
Qx = stl(x,15)$time.series[,1]
tsplot(Qx, main="Seasonal Component", ylab=" % rooms",
ylim=c(-4.7,4.7), col=gray(.7))
text(Qx, labels=1:4, col=c(3,4,2,6), cex=.8)
par(mfrow=c(2,1))
tsplot(cardox, col=4, ylab=expression(CO[2]))
title("Monthly Carbon Dioxide Readings - Mauna Loa Observatory",
cex.main=1)
tsplot(diff(diff(cardox,12)), col=4,
ylab=expression(nabla~nabla[12]~CO[2]))
acf2(diff(diff(cardox,12)))
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21]
ACF -0.3 -0.03 -0.06 0.00 0.00 -0.01 0.03 -0.01 0.05 -0.03 0.19 -0.47 0.08 0.08 0.00 0.02 -0.05 0.04 -0.04 -0.01 -0.02
PACF -0.3 -0.13 -0.12 -0.07 -0.04 -0.03 0.02 0.01 0.07 0.02 0.22 -0.39 -0.20 -0.03 -0.06 -0.02 -0.07 0.00 -0.01 -0.03 0.02
[,22] [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37] [,38] [,39] [,40] [,41]
ACF 0.02 -0.04 0.00 0.03 -0.06 0.04 -0.03 0.06 -0.03 -0.04 0.04 0.01 0.02 0.00 0.01 0.04 -0.01 -0.01 -0.01 -0.01
PACF 0.00 0.08 -0.25 -0.15 -0.08 -0.04 -0.04 -0.01 0.02 -0.06 -0.02 0.02 0.03 0.07 -0.17 -0.02 -0.06 -0.02 -0.06 -0.03
[,42] [,43] [,44] [,45] [,46] [,47] [,48]
ACF -0.03 0.03 0.03 -0.05 0 0.02 -0.02
PACF -0.04 -0.11 0.02 -0.05 0 0.07 -0.14
sarima(cardox, p=0,d=1,q=1, P=0,D=1,Q=1,S=12)
initial value -0.826338
iter 2 value -1.059073
iter 3 value -1.093845
iter 4 value -1.116555
iter 5 value -1.124382
iter 6 value -1.126345
iter 7 value -1.127354
iter 8 value -1.127953
iter 9 value -1.127984
iter 10 value -1.127985
iter 10 value -1.127985
iter 10 value -1.127985
final value -1.127985
converged
initial value -1.144615
iter 2 value -1.148048
iter 3 value -1.148645
iter 4 value -1.149895
iter 5 value -1.150013
iter 6 value -1.150021
iter 7 value -1.150021
iter 7 value -1.150021
iter 7 value -1.150021
final value -1.150021
converged
<><><><><><><><><><><><><><>
Coefficients:
Estimate SE t.value p.value
ma1 -0.3869 0.0377 -10.2624 0
sma1 -0.8655 0.0183 -47.2846 0
sigma^2 estimated as 0.0980908 on 766 degrees of freedom
AIC = 0.5456475 AICc = 0.545668 BIC = 0.5637873
sarima(cardox, 1,1,1, 0,1,1,12)
initial value -0.827261
iter 2 value -1.034332
iter 3 value -1.066295
iter 4 value -1.094823
iter 5 value -1.108013
iter 6 value -1.115246
iter 7 value -1.116173
iter 8 value -1.120248
iter 9 value -1.120892
iter 10 value -1.121657
iter 11 value -1.122186
iter 12 value -1.124590
iter 13 value -1.125269
iter 14 value -1.125654
iter 15 value -1.125685
iter 16 value -1.125685
iter 17 value -1.125687
iter 18 value -1.125689
iter 19 value -1.125689
iter 20 value -1.125689
iter 20 value -1.125689
iter 20 value -1.125689
final value -1.125689
converged
initial value -1.146682
iter 2 value -1.150731
iter 3 value -1.152123
iter 4 value -1.152815
iter 5 value -1.153157
iter 6 value -1.153220
iter 7 value -1.153266
iter 8 value -1.153337
iter 9 value -1.153352
iter 10 value -1.153359
iter 11 value -1.153384
iter 12 value -1.153388
iter 13 value -1.153390
iter 14 value -1.153390
iter 14 value -1.153390
iter 14 value -1.153390
final value -1.153390
converged
<><><><><><><><><><><><><><>
Coefficients:
Estimate SE t.value p.value
ar1 0.2203 0.0894 2.4660 0.0139
ma1 -0.5797 0.0753 -7.7030 0.0000
sma1 -0.8656 0.0182 -47.5947 0.0000
sigma^2 estimated as 0.09742764 on 765 degrees of freedom
AIC = 0.541514 AICc = 0.5415549 BIC = 0.5657004
sarima.for(cardox, 60, 1,1,1, 0,1,1,12)
$pred
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
2023 422.5365 423.2516 422.6841 420.8017 418.7844 417.4195 417.6341 419.1994 420.6362
2024 421.8604 422.7144 423.3343 424.7951 425.4936 424.9223 423.0392 421.0216 419.6567 419.8713 421.4367 422.8734
2025 424.0976 424.9517 425.5715 427.0323 427.7308 427.1595 425.2764 423.2588 421.8940 422.1085 423.6739 425.1106
2026 426.3348 427.1889 427.8088 429.2695 429.9680 429.3968 427.5136 425.4961 424.1312 424.3458 425.9111 427.3479
2027 428.5720 429.4261 430.0460 431.5067 432.2052 431.6340 429.7509 427.7333 426.3684 426.5830 428.1483 429.5851
2028 430.8092 431.6633 432.2832
$se
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
2023 0.3121340 0.3706892 0.4100237 0.4437896 0.4747345 0.5036938 0.5310578 0.5570754 0.5819301
2024 0.6057658 0.6286983 0.6508233 0.6839230 0.7112115 0.7366194 0.7609956 0.7845756 0.8074590 0.8297096 0.8513785 0.8725094
2025 0.8931404 0.9133056 0.9330350 0.9616380 0.9859772 1.0090191 1.0313946 1.0532622 1.0746779 1.0956735 1.1162740 1.1365011
2026 1.1563743 1.1759118 1.1951299 1.2221148 1.2455209 1.2678703 1.2896982 1.3111337 1.3322181 1.3529726 1.3734132 1.3935539
2027 1.4134077 1.4329864 1.4523012 1.4786701 1.5018653 1.5241385 1.5459682 1.5674673 1.5886697 1.6095916 1.6302447 1.6506393
2028 1.6707850 1.6906907 1.7103647
abline(v=2018.9, lty=6)
##-- for comparison, try the first model --##
sarima.for(cardox, 60, 0,1,1, 0,1,1,12) # not shown
$pred
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
2023 422.4940 423.1899 422.6182 420.7348 418.7172 417.3524 417.5670 419.1324 420.5691
2024 421.7933 422.6474 423.2671 424.7222 425.4182 424.8465 422.9630 420.9455 419.5807 419.7953 421.3606 422.7974
2025 424.0216 424.8757 425.4954 426.9505 427.6464 427.0747 425.1913 423.1738 421.8089 422.0235 423.5889 425.0257
2026 426.2498 427.1039 427.7237 429.1788 429.8747 429.3030 427.4196 425.4020 424.0372 424.2518 425.8172 427.2539
2027 428.4781 429.3322 429.9519 431.4070 432.1030 431.5312 429.6478 427.6303 426.2655 426.4800 428.0454 429.4822
2028 430.7063 431.5605 432.1802
$se
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
2023 0.3131945 0.3673801 0.4145425 0.4568619 0.4955806 0.5314860 0.5651148 0.5968518 0.6269843
2024 0.6557337 0.6832745 0.7097474 0.7473726 0.7784772 0.8083858 0.8372267 0.8651067 0.8921157 0.9183308 0.9438180 0.9686348
2025 0.9928314 1.0164523 1.0395365 1.0715498 1.0989071 1.1255998 1.1516739 1.1771707 1.2021268 1.2265752 1.2505458 1.2740654
2026 1.2971587 1.3198479 1.3421537 1.3722431 1.3984561 1.4241867 1.4494607 1.4743014 1.4987305 1.5227677 1.5464314 1.5697383
2027 1.5927042 1.6153437 1.6376702 1.6670522 1.6930078 1.7185714 1.7437603 1.7685905 1.7930768 1.8172333 1.8410728 1.8646075
2028 1.8878489 1.9108076 1.9334937
trend = time(cmort); temp = tempr - mean(tempr); temp2 = temp^2
fit = lm(cmort~trend + temp + temp2 + part, na.action=NULL)
acf2(resid(fit), 52) # implies AR2
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21]
ACF 0.34 0.44 0.28 0.28 0.16 0.12 0.07 0.01 0.03 -0.05 -0.02 0.00 -0.04 -0.02 0.01 -0.05 -0.01 -0.03 -0.06 -0.03 -0.03
PACF 0.34 0.36 0.08 0.06 -0.05 -0.05 -0.02 -0.05 0.02 -0.06 0.00 0.06 -0.02 0.00 0.04 -0.08 0.00 -0.01 -0.06 0.03 0.02
[,22] [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37] [,38] [,39] [,40] [,41]
ACF -0.05 0.00 -0.02 -0.01 -0.05 0.03 -0.11 -0.02 -0.10 -0.02 -0.09 -0.03 -0.10 -0.08 -0.10 -0.07 -0.07 -0.05 -0.05 -0.04
PACF -0.03 0.05 -0.01 0.00 -0.06 0.06 -0.11 0.00 -0.05 0.05 -0.04 0.04 -0.06 -0.06 -0.05 0.03 -0.02 0.02 -0.01 0.00
[,42] [,43] [,44] [,45] [,46] [,47] [,48] [,49] [,50] [,51] [,52]
ACF -0.03 -0.04 0.05 0.00 0.04 0.08 0.07 0.06 0.05 0.07 0.02
PACF 0.00 -0.01 0.08 -0.01 -0.01 0.08 -0.01 -0.01 -0.03 0.03 -0.04
sarima(cmort, 2,0,0, xreg=cbind(trend, temp, temp2, part) )
initial value 1.849900
iter 2 value 1.733730
iter 3 value 1.701063
iter 4 value 1.682463
iter 5 value 1.657377
iter 6 value 1.652444
iter 7 value 1.641726
iter 8 value 1.635302
iter 9 value 1.630848
iter 10 value 1.629286
iter 11 value 1.628731
iter 12 value 1.628646
iter 13 value 1.628634
iter 14 value 1.628633
iter 15 value 1.628632
iter 16 value 1.628628
iter 17 value 1.628627
iter 18 value 1.628627
iter 19 value 1.628626
iter 20 value 1.628625
iter 21 value 1.628622
iter 22 value 1.628618
iter 23 value 1.628614
iter 24 value 1.628612
iter 25 value 1.628611
iter 26 value 1.628610
iter 27 value 1.628610
iter 28 value 1.628609
iter 29 value 1.628609
iter 30 value 1.628608
iter 31 value 1.628608
iter 32 value 1.628608
iter 32 value 1.628608
iter 32 value 1.628608
final value 1.628608
converged
initial value 1.630401
iter 2 value 1.630393
iter 3 value 1.630382
iter 4 value 1.630381
iter 5 value 1.630375
iter 6 value 1.630370
iter 7 value 1.630362
iter 8 value 1.630354
iter 9 value 1.630349
iter 10 value 1.630347
iter 11 value 1.630347
iter 12 value 1.630347
iter 13 value 1.630347
iter 14 value 1.630346
iter 15 value 1.630346
iter 16 value 1.630346
iter 17 value 1.630346
iter 17 value 1.630346
iter 17 value 1.630346
final value 1.630346
converged
<><><><><><><><><><><><><><>
Coefficients:
Estimate SE t.value p.value
ar1 0.3848 0.0436 8.8329 0.0000
ar2 0.4326 0.0400 10.8062 0.0000
intercept 3075.1482 834.7278 3.6840 0.0003
trend -1.5165 0.4226 -3.5881 0.0004
temp -0.0190 0.0495 -0.3837 0.7014
temp2 0.0154 0.0020 7.6117 0.0000
part 0.1545 0.0272 5.6803 0.0000
sigma^2 estimated as 26.01476 on 501 degrees of freedom
AIC = 6.130066 AICc = 6.130507 BIC = 6.196687
library(zoo)
lag2.plot(Hare, Lynx, 5) # lead-lag relationship
pp = as.zoo(ts.intersect(Lynx, HareL1 = lag(Hare,-1)))
summary(reg <- lm(pp$Lynx~ pp$HareL1)) # results not displayed
Call:
lm(formula = pp$Lynx ~ pp$HareL1)
Residuals:
Min 1Q Median 3Q Max
-31.107 -11.511 -2.193 7.579 45.632
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 15.84712 2.75817 5.746 1.29e-07 ***
pp$HareL1 0.27265 0.04727 5.768 1.17e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 16.25 on 88 degrees of freedom
Multiple R-squared: 0.2744, Adjusted R-squared: 0.2661
F-statistic: 33.27 on 1 and 88 DF, p-value: 1.172e-07
acf2(resid(reg)) # in Figure 5.19
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20]
ACF 0.73 0.30 -0.10 -0.35 -0.48 -0.43 -0.21 0.03 0.21 0.28 0.27 0.09 -0.16 -0.32 -0.35 -0.26 -0.10 0.14 0.34 0.37
PACF 0.73 -0.52 -0.18 -0.05 -0.25 0.08 0.13 -0.08 0.04 0.03 0.00 -0.24 -0.10 0.07 -0.13 0.04 0.06 0.08 0.05 -0.10
( reg2 = sarima(pp$Lynx, 2,0,0, xreg=pp$HareL1 ))
initial value 2.758050
iter 2 value 2.608114
iter 3 value 2.338840
iter 4 value 2.248812
iter 5 value 2.046836
iter 6 value 2.045593
iter 7 value 2.045218
iter 8 value 2.045114
iter 9 value 2.045060
iter 10 value 2.045055
iter 11 value 2.045055
iter 12 value 2.045055
iter 13 value 2.045055
iter 13 value 2.045055
iter 13 value 2.045055
final value 2.045055
converged
initial value 2.056766
iter 2 value 2.056665
iter 3 value 2.056583
iter 4 value 2.056581
iter 5 value 2.056580
iter 6 value 2.056580
iter 6 value 2.056580
iter 6 value 2.056580
final value 2.056580
converged
<><><><><><><><><><><><><><>
Coefficients:
Estimate SE t.value p.value
ar1 1.3258 0.0732 18.1184 0.0000
ar2 -0.7143 0.0731 -9.7689 0.0000
intercept 25.1319 2.5469 9.8676 0.0000
xreg 0.0692 0.0318 2.1727 0.0326
sigma^2 estimated as 59.57091 on 86 degrees of freedom
AIC = 7.062148 AICc = 7.067377 BIC = 7.201026
$fit
Call:
arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, Q), period = S),
xreg = xreg, transform.pars = trans, fixed = fixed, optim.control = list(trace = trc,
REPORT = 1, reltol = tol))
Coefficients:
ar1 ar2 intercept xreg
1.3258 -0.7143 25.1319 0.0692
s.e. 0.0732 0.0731 2.5469 0.0318
sigma^2 estimated as 59.57: log likelihood = -312.8, aic = 635.59
$degrees_of_freedom
[1] 86
$ttable
Estimate SE t.value p.value
ar1 1.3258 0.0732 18.1184 0.0000
ar2 -0.7143 0.0731 -9.7689 0.0000
intercept 25.1319 2.5469 9.8676 0.0000
xreg 0.0692 0.0318 2.1727 0.0326
$ICs
AIC AICc BIC
7.062148 7.067377 7.201026
prd = Lynx - resid(reg2$fit) # prediction (resid = obs - pred)
prde = sqrt(reg2$fit$sigma2) # prediction error
tsplot(prd, lwd=2, col=rgb(0,0,.9,.5), ylim=c(-20,90), ylab="Lynx")
points(Lynx, pch=16, col=rgb(.8,.3,0))
x = time(Lynx)[-1]
xx = c(x, rev(x))
yy = c(prd - 2*prde, rev(prd + 2*prde))
polygon(xx, yy, border=8, col=rgb(.4, .5, .6, .15))
mtext(expression(""%*% 10^3), side=2, line=1.5, adj=.975)
legend("topright", legend=c("Predicted", "Observed"), lty=c(1,NA),
lwd=2, pch=c(NA,16), col=c(4,rgb(.8,.3,0)), cex=.9)
t = seq(0, 24, by=.01)
X = cos(2*pi*t*1/2) # 1 cycle every 2 hours
tsplot(t, X, xlab="Hours")
T = seq(1, length(t), by=250) # observed every 2.5 hrs
points(t[T], X[T], pch=19, col=4)
lines(t, cos(2*pi*t/10), col=4)
axis(1, at=t[T], labels=FALSE, lwd.ticks=2, col.ticks=2)
abline(v=t[T], col=rgb(1,0,0,.2), lty=2)
x1 = 2*cos(2*pi*1:100*6/100) + 3*sin(2*pi*1:100*6/100)
x2 = 4*cos(2*pi*1:100*10/100) + 5*sin(2*pi*1:100*10/100)
x3 = 6*cos(2*pi*1:100*40/100) + 7*sin(2*pi*1:100*40/100)
x = x1 + x2 + x3
par(mfrow=c(2,2))
tsplot(x1, ylim=c(-10,10), main=expression(omega==6/100~~~A^2==13))
tsplot(x2, ylim=c(-10,10), main=expression(omega==10/100~~~A^2==41))
tsplot(x3, ylim=c(-10,10), main=expression(omega==40/100~~~A^2==85))
tsplot(x, ylim=c(-16,16), main="sum")
P = Mod(fft(x)/sqrt(100))^2 # periodogram
sP = (4/100)*P # scaled peridogram
Fr = 0:99/100 # fundamental frequencies
tsplot(Fr, sP, type="o", xlab="frequency", ylab="scaled periodogram",
col=4, ylim=c(0,90))
abline(v=.5, lty=5)
abline(v=c(.1,.3,.7,.9), lty=1, col=gray(.9))
axis(side=1, at=seq(.1,.9,by=.2))
par(mfrow=c(3,2), mar=c(1.5,2,1,0)+1, mgp=c(1.6,.6,0))
for(i in 4:9){
mvspec(fmri1[,i], main=colnames(fmri1)[i], ylim=c(0,3), xlim=c(0,.2),
col=rgb(.05,.6,.75), lwd=2, type="o", pch=20)
abline(v=1/32, col="dodgerblue", lty=5) # stimulus frequency
}
par(mfrow=c(3,1))
arma.spec(main="White Noise", col=4)
arma.spec(ma=.5, main="Moving Average", col=4)
arma.spec(ar=c(1,-.9), main="Autoregression", col=4)
par(mfrow=c(3,1))
tsplot(soi, col=4, main="SOI")
tsplot(diff(soi), col=4, main="First Difference")
k = kernel("modified.daniell", 6) # MA weights
tsplot(kernapply(soi, k), col=4, main="Seasonal Moving Average")
##-- frequency responses --##
par(mfrow=c(2,1))
w = seq(0, .5, by=.01)
FRdiff = abs(1-exp(2i*pi*w))^2
tsplot(w, FRdiff, xlab="frequency", main="High Pass Filter")
u = cos(2*pi*w)+cos(4*pi*w)+cos(6*pi*w)+cos(8*pi*w)+cos(10*pi*w)
FRma = ((1 + cos(12*pi*w) + 2*u)/12)^2
tsplot(w, FRma, xlab="frequency", main="Low Pass Filter")
library(fGarch)
library(tseries)
library(arfima)
library(tsDyn)
library(ggplot2)
library(dplyr)
######################### Warning from 'xts' package ##########################
# #
# The dplyr lag() function breaks how base R's lag() function is supposed to #
# work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or #
# source() into this session won't work correctly. #
# #
# Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
# conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop #
# dplyr from breaking base R's lag() function. #
# #
# Code in packages is not affected. It's protected by R's namespace mechanism #
# Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning. #
# #
###############################################################################
Attaching package: ‘dplyr’
The following objects are masked from ‘package:xts’:
first, last
The following objects are masked from ‘package:stats’:
filter, lag
The following objects are masked from ‘package:base’:
intersect, setdiff, setequal, union
renv::init()
renv: Project Environments for R
Welcome to renv! It looks like this is your first time using renv.
This is a one-time message, briefly describing some of renv's functionality.
renv will write to files within the active project folder, including:
- A folder 'renv' in the project directory, and
- A lockfile called 'renv.lock' in the project directory.
In particular, projects using renv will normally use a private, per-project
R library, in which new packages will be installed. This project library is
isolated from other R libraries on your system.
In addition, renv will update files within your project directory, including:
- .gitignore
- .Rbuildignore
- .Rprofile
Finally, renv maintains a local cache of data on the filesystem, located at:
- "~/.cache/R/renv"
This path can be customized: please see the documentation in `?renv::paths`.
Please read the introduction vignette with `vignette("renv")` for more information.
You can browse the package documentation online at https://rstudio.github.io/renv/.
y
- "~/.cache/R/renv" has been created.
- Linking packages into the project library ... [73/128] [74/128] [75/128] [76/128] [77/128] [78/128] [79/128] [80/128] [81/128] [82/128] [83/128] [84/128] [85/128] [86/128] [87/128] [88/128] [89/128] [90/128] [91/128] [92/128] [93/128] [94/128] [95/128] [96/128] [97/128] [98/128] [99/128] [100/128] [101/128] [102/128] [103/128] [104/128] [105/128] [106/128] [107/128] [108/128] [109/128] [110/128] [111/128] [112/128] [113/128] [114/128] [115/128] [116/128] [117/128] [118/128] [119/128] [120/128] [121/128] [122/128] [123/128] [124/128] [125/128] [126/128] [127/128] [128/128] Done!
- Resolving missing dependencies ...
# Installing packages --------------------------------------------------------
The following package(s) will be updated in the lockfile:
# CRAN -----------------------------------------------------------------------
- abind [* -> 1.4-5]
- arfima [* -> 1.8-1]
- astsa [* -> 2.1]
- backports [* -> 1.4.1]
- base64enc [* -> 0.1-3]
- boot [* -> 1.3-28]
- brio [* -> 1.1.5]
- broom [* -> 1.0.5]
- bslib [* -> 0.7.0]
- cachem [* -> 1.0.8]
- callr [* -> 3.7.6]
- car [* -> 3.1-2]
- carData [* -> 3.0-5]
- cli [* -> 3.6.2]
- codetools [* -> 0.2-18]
- colorspace [* -> 2.1-0]
- cpp11 [* -> 0.4.7]
- crayon [* -> 1.5.2]
- curl [* -> 5.2.1]
- cvar [* -> 0.5]
- desc [* -> 1.4.3]
- deSolve [* -> 1.40]
- diffobj [* -> 0.3.5]
- digest [* -> 0.6.35]
- dplyr [* -> 1.1.4]
- dynlm [* -> 0.3-6]
- ellipsis [* -> 0.3.2]
- evaluate [* -> 0.23]
- fansi [* -> 1.0.6]
- farver [* -> 2.1.1]
- fastICA [* -> 1.2-4]
- fastmap [* -> 1.1.1]
- fBasics [* -> 4032.96]
- fGarch [* -> 4033.92]
- fontawesome [* -> 0.5.2]
- foreach [* -> 1.5.2]
- forecast [* -> 8.22.0]
- fracdiff [* -> 1.5-3]
- fs [* -> 1.6.4]
- gbutils [* -> 0.5]
- generics [* -> 0.1.3]
- ggplot2 [* -> 3.5.1]
- glue [* -> 1.7.0]
- gss [* -> 2.2-7]
- gtable [* -> 0.3.5]
- highr [* -> 0.10]
- htmltools [* -> 0.5.8.1]
- isoband [* -> 0.2.7]
- iterators [* -> 1.0.14]
- jquerylib [* -> 0.1.4]
- jsonlite [* -> 1.8.8]
- knitr [* -> 1.46]
- labeling [* -> 0.4.3]
- lattice [* -> 0.20-45]
- lifecycle [* -> 1.0.4]
- lme4 [* -> 1.1-35.3]
- lmtest [* -> 0.9-40]
- ltsa [* -> 1.4.6]
- magrittr [* -> 2.0.3]
- MASS [* -> 7.3-55]
- Matrix [* -> 1.6-5]
- MatrixModels [* -> 0.5-3]
- memoise [* -> 2.0.1]
- mgcv [* -> 1.8-39]
- mime [* -> 0.12]
- minqa [* -> 1.2.6]
- mnormt [* -> 2.1.1]
- munsell [* -> 0.5.1]
- nlme [* -> 3.1-155]
- nloptr [* -> 2.0.3]
- nnet [* -> 7.3-17]
- numDeriv [* -> 2016.8-1.1]
- pbkrtest [* -> 0.5.2]
- pillar [* -> 1.9.0]
- pkgbuild [* -> 1.4.4]
- pkgconfig [* -> 2.0.3]
- pkgload [* -> 1.3.4]
- praise [* -> 1.0.0]
- processx [* -> 3.8.4]
- ps [* -> 1.7.6]
- purrr [* -> 1.0.2]
- quadprog [* -> 1.5-8]
- quantmod [* -> 0.4.26]
- quantreg [* -> 5.97]
- R6 [* -> 2.5.1]
- rappdirs [* -> 0.3.3]
- rbibutils [* -> 2.2.16]
- RColorBrewer [* -> 1.1-3]
- Rcpp [* -> 1.0.12]
- RcppArmadillo [* -> 0.12.8.2.1]
- RcppEigen [* -> 0.3.4.0.0]
- Rdpack [* -> 2.6]
- rematch2 [* -> 2.1.2]
- renv [* -> 1.0.7]
- rlang [* -> 1.1.4]
- rmarkdown [* -> 2.26]
- rprojroot [* -> 2.0.4]
- sandwich [* -> 3.1-0]
- sass [* -> 0.4.9]
- scales [* -> 1.3.0]
- SparseM [* -> 1.81]
- spatial [* -> 7.3-15]
- stabledist [* -> 0.7-1]
- stringi [* -> 1.8.3]
- stringr [* -> 1.5.1]
- strucchange [* -> 1.5-3]
- survival [* -> 3.2-13]
- testthat [* -> 3.2.1.1]
- tibble [* -> 3.2.1]
- tidyr [* -> 1.3.1]
- tidyselect [* -> 1.2.1]
- timeDate [* -> 4032.109]
- timeSeries [* -> 4032.109]
- tinytex [* -> 0.50]
- tsDyn [* -> 11.0.4.1]
- tseries [* -> 0.10-56]
- tseriesChaos [* -> 0.1-13.1]
- TTR [* -> 0.24.4]
- urca [* -> 1.3-4]
- utf8 [* -> 1.2.4]
- vars [* -> 1.6-1]
- vctrs [* -> 0.6.5]
- viridisLite [* -> 0.4.2]
- waldo [* -> 0.5.2]
- withr [* -> 3.0.0]
- xfun [* -> 0.43]
- xts [* -> 0.14.0]
- yaml [* -> 2.3.8]
- zoo [* -> 1.8-12]
The version of R recorded in the lockfile will be updated:
- R [* -> 4.1.2]
- Lockfile written to "~/water-quality-time-series/renv.lock".